#include "../includes/need-a-name.hh"
#include "includes/initcollecter-euler.hh" 
#include <string>

InitCollecterEuler::InitCollecterEuler(){
  set_parameter("cen1",0.0);
  set_parameter("cen2",0.0);
  set_parameter("cen3",0.0);
  set_parameter("slope1",0.0);
  set_parameter("slope2",0.0);
  set_parameter("slope3",0.0);  
  set_parameter("radius",0.0);
  
  set_parameter("rhoI",0.0);
  set_parameter("rhoII",0.0);    
  set_parameter("mom1I",0.0);
  set_parameter("mom1II",0.0);

  set_parameter("eI",0.0);
  set_parameter("eII",0.0);
  set_parameter("gamma",5.0/3.0);  
}
InitCollecterEuler::~InitCollecterEuler(){
}

void InitCollecterEuler::initial_function(double* location, double pillarlen, int currentdim, int index){
  double componenttmp;  
  double basisInside[memory::dgorder*memory::intergnumber];
  double locationtmp[gridsetting::dimension];
  if(memory::dgflag==0){
    for(int tmpi=0; tmpi<memory::phynumber;tmpi++){
      if(para["sod_shocktube"]==1.0)
	memory::variable[tmpi][index]=sod_shocktube(location, tmpi);
      if(para["blastwave"]==1.0)
	memory::variable[tmpi][index]=blastwave(location, tmpi);
      if(para["two_shockwave_1d"]==1.0)
	memory::variable[tmpi][index]=two_shockwave_1d(location, tmpi);      
    }
  }

  if(memory::dgflag==1){
    basis.getBasis(basisInside);
    for(int tmpl=0; tmpl<gridsetting::dimension; tmpl++)
      locationtmp[tmpl]=location[tmpl];
    for(int tmpi=0; tmpi<memory::phynumber;tmpi++)
      for(int tmpj=0; tmpj<memory::dgorder;tmpj++){
    	componenttmp=0.0;
    	for(int tmpk=0; tmpk<memory::intergnumber; tmpk++){
    	  locationtmp[currentdim]=location[currentdim]+(-0.5+tmpk/(memory::intergnumber-1.0))*pillarlen;
	  if(para["two_shockwave_1d"]==1.0)
	    componenttmp+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*two_shockwave_1d(locationtmp,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];	  
	  if(para["sod_shocktube"]==1.0)
	    componenttmp+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*sod_shocktube(locationtmp,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];
	  if(para["blastwave"]==1.0)
	    componenttmp+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*blastwave(locationtmp,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];
	}
	memory::variable[tmpj*memory::phynumber+tmpi][index]=componenttmp;
      }
    for(int tmpi=0; tmpi<memory::phynumber;tmpi++){
      if(para["two_shockwave_1d"]==1.0)
    	memory::variable[memory::dgorder*memory::phynumber+tmpi][index]=two_shockwave_1d(location,tmpi);
      if(para["sod_shocktube"]==1.0)
    	memory::variable[memory::dgorder*memory::phynumber+tmpi][index]=sod_shocktube(location,tmpi);      
      if(para["blastwave"]==1.0)
    	memory::variable[memory::dgorder*memory::phynumber+tmpi][index]=blastwave(location,tmpi);
    }
  }  
}

// Figure 14.4 of book: Leveque R. J.
double InitCollecterEuler::two_shockwave_1d(double* location, int phyi){
  double tmp=0.0;
  para["cen1"]=0.3;  
  para["rhoI"]=1.0;
  para["rhoII"]=2.0;    
  para["mom1I"]=1.0*3.0;
  para["mom1II"]=2.0*1.0; 
  para["pI"]=1.0;
  para["pII"]=1.0;
  if(phyi==0)
    tmp = (location[0]<para["cen1"]) ? para["rhoI"] : para["rhoII"];
  if(phyi==1)
    tmp = (location[0]<para["cen1"]) ? para["mom1I"] : para["mom1II"];
  if(phyi==2)
    tmp = (location[0]<para["cen1"]) ? std::pow(para["momI"],2)/para["rhoI"]/2.0+para["pI"]/(para["gamma"]-1.0) : std::pow(para["momII"],2)/para["rhoII"]/2.0+para["pII"]/(para["gamma"]-1.0);
  return tmp;
}

// 1-dimension: Figure 14.3 of book: Leveque R. J.
double InitCollecterEuler::sod_shocktube(double* location, int phyi){
  double cen[gridsetting::dimension];
  double slope[gridsetting::dimension];  
  double dotwithnormal=0;
  double tmp=0.0;
  para["cen1"]=0.5;
  para["cen2"]=0.5;
  para["cen3"]=0.5;
  para["slope1"]=1.0;
  para["slope2"]=1.0;
  para["slope3"]=1.0;    
  para["rhoI"]=3.0;
  para["rhoII"]=1.0;    
  para["pI"]=3.0;
  para["pII"]=1.0;
  if(gridsetting::dimension>=1){
    cen[0]=para["cen1"];
    slope[0]=para["slope1"];
  }
  if(gridsetting::dimension>=2){
    cen[1]=para["cen2"];
    slope[1]=para["slope2"];
  }
  if(gridsetting::dimension==3){
    cen[2]=para["cen3"];
    slope[2]=para["slope3"];    
  }
  // memory::variable[0][index] = (location[0]<para["cen1"]) ? para["rhoI"] : para["rhoII"];
  // memory::variable[1][index] = 0.0;
  // memory::variable[2][index] = (location[0]<para["cen1"]) ? std::pow(para["momI"],2)/para["rhoI"]/2.0+para["pI"]/(para["gamma"]-1.0) : std::pow(para["momII"],2)/para["rhoII"]/2.0+para["pII"]/(para["gamma"]-1.0);
  for(int tmpi=0; tmpi<gridsetting::dimension;tmpi++)
    dotwithnormal=dotwithnormal+slope[tmpi]*(location[tmpi]-cen[tmpi]);
  if(phyi==0)
    tmp = dotwithnormal>0? para["rhoI"] : para["rhoII"];
  //for(int tmpi=0; tmpi<gridsetting::dimension; tmpi++)
  //  memory::variable[tmpi+1][index] = 0.0;
  if(phyi>=1 && phyi<=gridsetting::dimension)
    tmp = 0.0;  
  if(phyi==gridsetting::dimension+1)
    tmp = dotwithnormal>0? para["pI"]/(para["gamma"]-1.0) : para["pII"]/(para["gamma"]-1.0);  
  
  return tmp;
}

double InitCollecterEuler::blastwave(double* location, int phyi){
  para["rhoI"]=1.0;
  para["pI"]=1.5;
  para["pII"]=1.0;
  para["radius"]=0.25;
  double distance=0;
  double tmp;
  if(phyi==0)
    tmp = para["rhoI"];
  if(phyi>=1 && phyi<=gridsetting::dimension)
    tmp = 0.0;
  for(int tmpi=0; tmpi<gridsetting::dimension; tmpi++){
    distance+=std::pow(location[tmpi]-0.5,2);
  }
  distance=std::sqrt(distance);  
  if(phyi==gridsetting::dimension+1)
    tmp = distance<=para["radius"] ? para["pI"]/(para["gamma"]-1.0) : para["pII"]/(para["gamma"]-1.0);
  return tmp;
}

